#!/usr/bin/env python3 
# -*- coding: utf-8 -*-
#-------------------------------------------------------
#	FileName	: spec.py
#	Author		：hpy
#	Date		：2024年03月02日
#	Description	：
#-------------------------------------------------------
from sys import argv
import numpy as np
import math 
from scipy.signal import find_peaks

#-------------------------------------------------------
#	从txt文件加载光谱数据
#   文件格式：第一列x  第二列y
#-------------------------------------------------------
def loadreadtxt(file_name):
    data = []
    file = open(file_name,'r')  #打开文件
    file_data = file.readlines() #读取所有行
    for row in file_data:
        row = row.strip("\n")
        row = row.strip("\r")
        row = row.strip()
        tmp_list = row.split(' ') #按‘，’切分每行的数据
        # tmp_list[-1] = tmp_list[-1].replace('\n','') #去掉换行符 \n
        # print(tmp_list)
        data.append(tmp_list) #将每行数据插入data中

    return data 


# 矩阵遍历打印到终端
def mat_traverse(mat):
    rows, cols = mat.shape              # 获取矩阵的行数和列数
    for i in range(rows):               # 按照行来遍历
        for j in range(cols):           # 对第i行进行遍历
            print('%.3f'% mat[i, j], end=' ')
            # print(mat[i, j], end=' ')   # 输出第i行第j列元素（读）
        print("")


# def find_nearest(a, a0):
#     "Element in nd array `a` closest to the scalar value `a0`"
#     idx = np.abs(a - a0).argmin()
#     return a.flat[idx]


def find_nearest(a, a0):
    "Element in nd array `a` closest to the scalar value `a0`"
    b= np.flipud(a)  # a的逆序
    re = []
    for i in a0:
        idx = np.abs(a - i).argmin()
        r_idx = np.abs(b - i).argmin()
        ave = (a.flat[idx] + b.flat[r_idx])/2.0
        re.append(ave)
    return re


# PS 片峰位置
peaks_dst = [
    842.32 , 
    906.84 , 
    1028.43,
    1069.21 , 
    1154.55 , 
    1538.21 ,
    1601.31,
    2849.93,
    3001.84,
    3026.15,
    3060.21,
    3082.25
]


if __name__=="__main__":
    # 加载背景光谱
    back = np.mat( loadreadtxt('src/background.back') ,  dtype=np.float64 )
    spec = back[: , 0:1]
    #将目标数据转为矩阵 
    peaks_out = np.mat(peaks_dst).T
    # 获取光谱
    first = True 
    for arg in argv:
        if first:
            first = False
            continue

        
        # print (arg)
        # 加载能量数据
        data = np.mat( loadreadtxt(arg) ,  dtype=np.float64 )
        # 计算吸光度，使用吸光度来寻峰
        # data[: , 1:2] =  np.log(back[: , 1:2]/data[: , 1:2]) / np.log(10)
        data[: , 1:2] =  np.log(back[: , 1:2]/data[: , 1:2]) / np.log(10)

        # 矩阵转数组
        arr_spec = np.array(data[: , 1:2])
        # 将二维数组转一维数组
        arr_spec = arr_spec.reshape(-1)
        # print(arr_spec)
        # 寻峰的位置
        peaks, _ = find_peaks(arr_spec, distance=6)
        # 得到pk的值，波数位置
        pk = find_nearest( np.array(data[peaks , 0]), peaks_dst )
        # 将寻找的波数位置按列连接到矩阵
        peaks_out = np.c_[peaks_out , pk]
        # print(peaks_out) 

        # print(peaks)
        # for peaks_idx in peaks:
        #         print(data[peaks_idx , 0] ,data[peaks_idx , 1]  )

        # print(data[peaks , 0] )
        # print(find_nearest( np.array(data[peaks , 0]).reshape(-1), peaks_dst ) )
        # print(find_nearest( np.array(data[peaks , 0]), peaks_dst ) )
        # spec = np.c_[spec,data[: , 1:2]]

    # 打印结果
    # mat_traverse(peaks_out) 

    dlt_cm = peaks_out
    dlt_cm [:,1:] = dlt_cm [:,1:]  - np.mat(peaks_dst).T
    dlt_ave = np.average(dlt_cm [:,1:] ,axis=1 )
    
    # 将最后一列的平均值计算出来
    print("# cm-1","|find cm","|dlt wavenumber ave")
    mat_traverse(np.c_[peaks_out,dlt_ave]) 
    # print (dlt_ave)
    # 矩阵按行求方差
    # spec[: , 1:] = np.std(spec[: , 1:],axis=1)

    # 按行求平均值
    # ave = np.average(spec[: , 1:], axis=1) 
    # spec = np.c_[spec,ave]
    # print(ave)

    # mat_traverse(spec)

    # back = np.mat( loadreadtxt('src/background.back') ,  dtype=np.float64 )
    # data = np.mat( loadreadtxt(srcfile) ,  dtype=np.float64 )
    # ans  = np.copy(back)

    # 计算透过率
    # ans[: , 1:2] =data[: , 1:2] / back[: , 1:2] * 100 

    # 计算吸光度
    # ans[: , 1:2] = np.log(back[: , 1:2] / data[: , 1:2] ) / np.log(10)

    # 设置打印所有的矩阵数据
    # np.set_printoptions(threshold=np.inf)
    # print(spec)